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A fully nonlinear numerical simulation of two-dimensional Faraday waves between two 
incompressible and immiscible fluids is performed by adopting the phase-field method 
with the Cahn-Hilliard equation due to Jacqmin (J. Comput. Phys., vol. 155, 1999, 
pp. 96-127). Its validation is checked against the linear theory. In the nonlinear regime, 
qualitative comparison is made with an earlier vortex-sheet simulation of two dimensional 
Faraday waves by Wright , Yon & Pozrikidis (J. Fluid Mech., vol. 400, 2000, pp. 1-32). The 
vorticity outside the interface region is studied in this comparison. The period tripling 
state, which is observed in the quasi-two dimensional experiment by Jiang, Perlin & 
Schultz (J. Fluid Mech., vol. 369, 1998, pp. 273-299), is successfully simulated with the 
present phase-field method. 
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1. Introduction 

Faraday waves, which typically refer to complex patterns of standing waves on a fluid 
surf ace in an oscil lating container, are among the classical problems of fluid mechan- 
ics ( Faradavl 11831 ). The phenomenon has been a representative example of paramet- 
ric instabilities (jMiles fe Henderson! 1990). As is often the case with fluid mechanics, 
even classical phenomena are often not well understood in nonlinear regimes. Indeed, 
continuing experiments on Faraday waves beyond the linear regimes reveal intrigu- 
ing new features, which include snake-li ke structures in drop-confined Faraday waves 
(IPucci. Fort. Ben Armar fe Couderl 120111) . a turbulent state m ediated by defects of the 
pattern (jShani. Cohen fe Finebereil2010l ) and so-called oscillons (|Arbell fc Fineberd2000l) . 

These surprising findings would probably be outsid e of the applicable ra nge of the 
weakly nonlinear theories, such as the one developed bv lChen fc Vinalsl (|1999i) . on selec- 
tion of various patterns of Faraday waves. Hence, to understand these phenomena, fully 

nonlinear numerical simulations of the Faraday systems, which c an be complementary to 

labora tory experiments, play an indispensable role as discussed in lKitvk. Embs. Mekhonoshin fc Wagner 

( 20051 ). Perhaps the first such simulation, in which the motions of both the top and bot- 

torn fl uids are simultaneously simulated, was performed recently bv lPerinet. Juric fc Tuckerman 
(|2009l ). This motivates our present study. 

In these fully nonlinear simulations, the inevitable difficulties are how to represent 
numerically the interface between the two immiscible flui ds, how to follow its motion 
and h ow to calculate its influence on the bulk fluids. In IPerinet. Juric fc Tuckerman 
( 20091 ) . the interface is modelled by triangular elements whose apices are advected ver- 
tically and data on the interface elements are copied to the Eulerian grids in the bulk 
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of the fluids according to the recipe of the immersed boundary method (|Peskmlll977l ). 
Since numerical modelling of the interface involves various assumptions, its validation 
is required against laboratory experimental data or other data independent of the mod- 
elling. In the linear regime, the an alytical result of the linear theory of Faraday waves 
due to Kumar fe Tuckermanl (J1994') provides reliable data for comparison, as used in 
iPerinet. Juric fc Tuckermanl (2009). In nonlinear reg imes, comparison with experimenta l 
data is crucial. Remarkably, the simulation result bv lPerinet. Juric fc Tuckermanl ( 2009T ) 
is in perfect agreement with the experimental data in the nonlinear regime conducted 
bv lKitvk. Embs. Mekhonoshin fc Wagnerl ( 20051 ). where the top fluid's motion cannot be 
neglected. 

The further task of such validated numerical simulations is to investigate data not 
easily accessible in experiments, such as nonlinear energy transfers between modes. For 
this purpose, in our opinion, it is important to have at least two validated numerical 
simulations with independent interface modelling and make sure that these simulations 
give consistent results. 

The aim of this study is to develop another nonlinear si mulation method of Faraday 
waves with an alternative interface model to that devised bv lPerinet. Juric fc Tuckermanl 
( 20091 ). We here apply to the Faraday wa ve problem the phase-field modelling with the 
Cahn-Hilliard equation for binary fluids ( Jacqminlll999l ). The phase- field method, like 
front-tracking methods, volume-of-fluid methods and level-set methods, easily allows sit- 
uations where th e interface becomes a multivalued function of the horizont al coordinates. 
See, for example. ICelani. Mazzino. Muratore-Ginanneschi fc Vozellal ( 20091 ) for recent ap- 
plication of the phase-field method to Ray leigh- Taylor instability and to other flows in 
the references therein. To our knowledge, the method hat not been applied to the Faraday 
wave problem. 

In this paper we focus on Faraday waves in two spatial dimensions (2D) to explore 
the capabilities of the phase-fiel d method in a simp l e setting. Earlier num e rical stud- 
ies of 2D F araday waves include IChen fc Wul (|2000f ) . iMurakami fc Chikanol (J200lh and 
Ubal (2003), where the fluid dynamical equations for the bottom fluid are solved but the 
top fluid's motion is neglected in contrast to the present simulation. Another numeri- 
cal approac h using different formalisms (the boundary integral and the vortex sheet) is 
explored by I Wright, Yon fc Pozrikidisl (J2000l) . 

Here we proceed as follows. In the linear r egime, we compare quantitat ively the phase- 
field simulation with the linear t h eory by iKumar fc Tuckermanl ( 19941 ) for validation 



Perinet. Juric fc Tuckermanl ( 20091) . In the nonlinear regime, unfortunately suit 



as in 

able experimental data in 2D for comparison with our simulation are not available. 
However, we qualitatively make com parison with the simulation of 2 D Faraday waves 
with the vortex-sheet formulation bv I Wright. Yon fc Pozrikidisl ( 20001 ) in the regime of 
plume formation, where the interface becomes a multivalued function. Finally, we present 
simulation results concerning the period tripling state, where the oscillation period of 
the Faraday waves becomes three t i mes t he basic period. This state is found experi- 
mentally by Jiang. Perlin fc Schultzl ( 19981) in quasi- two dimensional Faraday waves, by 
iDas fc Hopfingerl (120081) in a xisymmetri c thre e-dimensional Faraday waves, and numer- 



ically also by Wrig ht. Yon fc Pozrikidisl (2000). Its underlying physics is not yet under- 
stood. 

The organization of the paper is the following. In §2, we describe the basic fluid dy- 
namical the equations and our phase-field model of the Faraday waves. Our numerical 
method for solving the equations is explained in §2. We then make a quantitative com- 
parison of its simulation results with the linear analysis in §3. In §4, we present results of 
the phase-field simulation in nonlinear regimes. The first simulation shows plume forma- 
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Figure 1. Faraday wave setting in two-dimensional space. 

tion, in which the interface overturns. The second one concerns the period tripling state. 
Concluding remarks are made in the last §5. 



2. Equations and numerical method 

2.1. Equations 

We consider Faraday waves between two immiscible fluids in two spatial dimensions. The 
fluid dynamical equations are the incompressible Navier-Stokes equations 



P 



9u , _. 



-Vp + pG + V ■ i](Vu + Vu T ) + s, 



V -u = 0, 



(2.1) 
(2.2) 



where p is the pressure and u is the velocity. Here s is the surface tension force, which 
will be discussed below. The density p and viscosity 77 are distinct constants for each 
fluid, which are denoted as pt,Pb and r?t,?7b respectively, for the top and bottom fluids 
{pt < Pb)- Finally, G is the gravitational acceleration, incorporating the vibration forcing 
with amplitude a and angular frequency lo, which is written as 



G = (—g + a cos wi)e 2 



(2.3) 



Here in two dimensions, e z — (0, 1). The temporal period of the vibration forcing is 
denoted by T v = 2tt/lo, which will be used often in what follows. 

The equations are solved with the boundary conditions of the Faraday setting as de- 
picted in figure [1] In the horizontal direction (coordinate x), we assume periodic bound- 
ary conditions with length L x . For the vertical direction (coordinate z), no-slip boundary 
conditions at the boundaries z — 0,L Z are assumed, that is, u(x, 0, t) = u(x,L z ,t) = 0. 
The interface position between the top and bottom fluids, denoted as z = ((x,t), obeys 
the kinematic boundary condition. In terms of this £(x, £), the density p and viscosity 7/ 
in (12.11) are written as function of z: 



(p. v) 



(p t ,Vt) z>C(x,t), 
(Pb,rjb) z^((x,t). 



(2.4) 
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In this sharp interface formulation, the density and viscosity vary discontinuously as a 
function of z, which makes numerical simulations in this form very difficult in practice. 
One standard way to overcome this in numerical simulations is to model the sharp inter- 
face by a diffuse one (| Anderson. McFadden fe Wheelerlll998l ). The phase-field modelling 
is one such diffuse-interface method, which we will ex plain below. 

We adapt the phase-field modelling of binary fluids ( Jacqmmlll999l ) to Faraday waves. 
The top and bottom fluids are here indicated by the phase variable values <f> = I and <j) = 
— I, respectively. The interface is modelled as the region where (j> changes continuously 
from —I to I. The phase evolves to minimize the free energy functional of the phase <f>: 



F[ct>(x,t)} = ±J 



i(^ 2 -l) 2 + l^ 



dx. 



(2.5) 



Accordingly, the equation of the phase variable is the Cahn-Hilliard equation with the 
advection term 



SF A 



(2.6) 

(2.7) 



where \x is the chemical potential, A is the magnitude of the free energy and 7 is the mo- 
bility. We assume that the mobility 7 is constant. The parameter e controls the thickness 
of the smoothed interface. These 7 and e are adjustable parameters of the phase-field 
modelling. How to determine them is discussed at the end of section [3] In terms of (j), 
the density and viscosity are expressed as 



Pt + Pb . Pt - Pb 



V 



Vt + Vb . Vt - Vb 



(2. 



In the framework of the phase-field modelling, the surface tension force written as s 
in (|2.I[) is given by 

s = fj,V(j) (2.9) 

( Jacaminll999tlCelani. Mazzino. Muratore-Ginanneschi fc Vozeiral2009l ) . The surface ten- 



sion a used in the sharp interface model is related to the phase-field simulation parameters 

as 



2V2A 



(2.10) 



for planar interfaces. In the following we assume that the correspondence (|2.10|) is valid 
for non-planar interfaces. We use the following boundary conditions of the phase variable 
on the top and bottom walls z — Q,L Z : 

e 2 -V0 = O, e z -V/i = (2.11) 

( Jacqminl 11999). In the horizontal direction, the boundary condition is peri odic. For 



furthe r details of the derivation of the phase-field model, we refer the reader to Uacamin 
(Il999h . 



To summarise this subsection, the equations to be solved numerically are the incom- 
pressible Navier-Stokes equations (|2.ip and (12.21) and the Cahn-Hilliard equation (12.6[) . 
The boundary conditions on the velocity are periodic in the horizontal direction and no- 
slip on the top and bottom walls in the two-dimensional space. For the phase variable, 
the boundary conditions are periodic in the x-direction and (|2.11|) in the z-direction. 
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2.2. Numerical method 

We start with the discretization of the C ahn-Hilliard equation. The computa tional mesh is 
a standard staggered arrangement as in lPerinet. Juric fc Tuckermanl (|2009l ) . We Fourier- 
expand the phase variable <j> in the periodic ^-direction as 



(x,z,t) 






3 



4>(k,z,t) e lfc ' 



(2.12) 



The wavenumber is given by kj = (2ir/L x )j (j ~ —N x /2 + 1, ...,N x /2) , where N x is 
the number of grid points in the x-coordinate. Then the Cahn-Hilliard equation in the 
(k, z) space becomes 



i:)t 



+ (u-V)(j) =jA- 



dz- 



i 



dz 



(J>-<P 3 ) 



(2.13) 



where the cubic term </> 3 is the Fourier transform of the product </> 3 calculated in the 
physical space (x, z). No aliasing error is removed in this calculation of 4> 3 . 

The advection term (u • V)(/> is calculated as follows. First it is calculated in the 
physical space by using an elementary discretization in the staggered mesh as 



[(«-V)fl, m = 



/; " ^dx 

Ul-\/2,m + u l+l/2. 



k+1, m — <t>l-l, m w l, m-1/2 + w l, m+1/2 4>l, m+1 ~ 4>l, ro-1 



2Ax 



2A; 



(2.14) 



In this notation u t _ 



1/2, 



denotes the velocity on the cell face (xz_i/2,Zm) and 



denote the phase variable on the cell centre {xi, z m ). The indices run as I = 1, • • • ,N X 
and m = !,••■ ,N Z (N z is the number of grid points on the ^-coordinates). In the 
denominator Aa; and Az are the grid spacings of the x and z coordinates. The Fourier 
transform of (|2.14p gives (u • V)0. Again no aliasing error is removed. 

Equation (12.1311 is discretized in time by using the semi-implicit stabilized scheme due 
to lEvrel(|l997t ) as 



(n+1) 



At 



7A ■ 



(«(") • V)</>( r 

,2 x2 



dz 



/>(»+!) 



-k 2 



dz 2 



(0 n) - 20 (n+1) - (^T)A 



(2.15) 



Here (f>( l > denotes the data at the i-th time step. On the right hand side, the fourth- 
order derivative term is treated fully implicitly. The second-order term is handled semi- 
implicitly, w hich i s split into two terms: one involves 30 ( ™- ) and the other 2^ n+1 ' as 
proposed by [Eyre] (|19971 ). The z-derivatives in (|2.15[) are evaluated as 



n 

= -7—4 ( 4>j,m+2 ~ ^4>j,m+l + &4>j,m ~ &<j>j,m-l 



a? 



(2.16) 
y„. 2 I - (2-17) 



6 

where 
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<j t m denotes cf>(kj,mAz). We use the bi-conjugate gradient stabilized (BiCGSTAB) 
method (|Saadl ll996) to solve the equations (|2.15j) . The boundary conditions of (j> in the 
z-direction are given in (|2.11[) . 

Now we turn to discretization of the incompressible Navicr-Sto kes equations (12.11) and 
( 2.21) . We follow mostly the discretization method proposed by Perinet. Juric fc Tuckerman 
( 20091 ) except for the surface tension force. In order to describe how we calculate the 



surface tension force in t he present phase-field modell i ng, w e repeat here some of the 
equations given in §3.3 of IPerinet. Juric fc Tuckerman! (2009). The method starts with 
the intermediate velocity u from u^ n ' at the n-th time step, 



,.(«) 



At 



-(«(") • V)u (n) 



1 



n(«+l) 



-V-»j< n+1 )(VM + V 



u T ) 



(2.18) 



where p^ and t/™) are here determined by the phase <f>( n > as ()2.8[) . Yet another interme- 
diate velocity u* involves the gravity and the vibration forces G and the surface tension 
force s: 



u — u 



At 

Since this s is represented as 
nent first as 



= G (n+1) 



s 



(ra+1) 



1 



Jn+1) 



Jn+1) 



TyVp 



(n) 



(2.19) 



in the phase-field method, we discretize each compo- 



(sJ n) 



An) 
In) Vl+l.m 



b (n) 



Jl,m r~L,m 2 Ax 

Then, the vectors s on the staggered grids are 



(s ) {n) 



An) 

(n) yj,m+l 



J l.rn-1 



2Al 



(8x) 



(n) 



\i<x)i + i^ m f \ b x)i^ m 



(«*) 



(n) 



\(n) 



(*.)SUi + («.), 1 



(n) 



Here the chemical potential (|2.7j) is discretized as 



(2.20) 



(2.21) 



J l+l,m 



2^ n) 

rirn 



+ 



b (n) 



An) 



20, (n) 



l.m- 



MS = A 



Finally, the velocity at the (n + l)-th time step is obtained by 
u ("+i) _ M * i 



& 



b (n) 
l.m 



At 



Jn+1) 



T)V[P 



,("+!) 



P 



(a) 



where the pressure p(" +1 ) is determined, by demanding V • u^ n+1 ^ = 0, as 
V -u* 



At 



p(n+l) \f 



P 



(n) 



(2.22) 



(2.23) 



(2.24) 



Again the BiCGSTAB method is employed to solve (J2.18I) and (I2.24[). The boundary con- 
dition s for the velocity and the intermediate one u are the same as in lPerinet. Juric fc Tuckerman 

( 20091 ). namely, u(x, z — 0oiL z ,t) = 0, e z • u^ n+Vl {x,z — or L ? ,t) = and e z • 

Vj^" + 1 )(x, z = 0orL z ,t)— p( n \x, z — 0orL 2 ,£)] = 0. As proposed bv lPerinet. Juric fc Tuckerman 
(|2009f ). in solving the equations (12.181) . we use the latest updated velocity component of 
u to compute the other component. In this process, the order (which component of u is 
computed first and last) is also interchanged at each time step to ensure symmetry. 
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Figure 2. Critical c urve of the vibration amplitude from the linear theory 

(Ku mar fc Tuckermanlll994l ). Points are the critical values determined with the phase- field 

simulation. 



3. Comparison in the linear regime 

We compare our simulation result with the linear theory of Kumar & Tuckcrman 
1994). Wc follow here again the validation method proposed bv lPerinet. Juric fc Tuckermanl 



2009). 



In the linear theory elaborated by Kumar fc Tuckermanl ( 1994I ). the Faraday wave 
problem is formulated as a finite-depth, viscous binary fluid system with a sharp interface 
representation. This provides the critical value of the vibration forcing amplitude a c 
as eigenvalues for a given perturbation e ikx . The Floquet modes f q (k) of the interface 
position z — C( x, t) can be calculated as eige nvectors. In other words, the spatio-temporal 
variation readsdKumar fc Tuckermanlll994h 



C(x,t)(xe ikx e (/5+iQw) * 



JV F 

£ 

q=-N F 



/«(*) 



^iqcjt 



(3.1) 



Recall that to is the angular frequency of the vibration forcing (|2.3p . Here the num- 
ber of Floquet modes JVp is infinite in theory but it is finite for numerical calcula- 
tions. The exponent (3 + iauj is the Floquet exponent. For the critical perturbations /3 
is zero. The harmonic a = and s ubharmonic a = 1/2 cases are considered here as in 
Perinet. Juric fc Tuckermanl (|2009l ). 

We numerically solve this eig envalue problem (for the precise form of the matrices, 



sec 
at 



Kumar & Tuckerman (1994)). The Mathematica script that we use here is available 



http://www.kyoryu.scphys.kyoto-u . ac .jp/~ltakeshi/kt94. The parameters here are the 
same as in IPerinet. Juric fc Tuckermanl ( 2009) in their validation with the linear theory. 
The parameter values are N F = 10, p h = 5.19933 x 10 2 kg m~ 3 , p t = 4.15667 x 10 2 kg 
m -3 j Vh = 3.908 x 10~ 5 Pa s, r/ t = 3.124 x 10~ 5 Pa s, a = 2.181 x 10~ 6 N m" 1 , g = 9.8066 
m s -2 , L z — 2.31 x 10 -4 m and u — 2tt x 10 2 s _1 . The unperturbed interface is in the 



centre of the container. The critical vibration forcing amplitude a c from the linear theory 
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k (mm 


) La (mm) 


a c /g (theory) 


a c /g (simu.) 


error (%) 


28.0 


0.224 


4.37 




4.65 


6.0 


48.0 


0.131 


12.5 




12.5 


0.0 


73.0 


0.0861 


28.5 




28.4 


0.49 


94.0 


0.0668 


51.0 




50.0 


2.0 



Table 1. Critical amplitude values a c of the linear theory calculation and the phase-field sim- 
ulation. Here the horizontal length of the simulation L x is taken to be equal to the wave- 
length of the perturbation 2-n/k for every case. The number of grid points of the simulation is 

N x x N z = 128 x 128 



as a function of the perturbation wavenumber k is plotted as a curve in figure [2j The 
Floquet modes f q are used later in figure |31 

With the phase-field simulation, critical values a c for four different fc's are determined, 
which are denoted as points in figure [2] Refer also to table Q] for a precise comparison 
with the values from the linear theory calculation. We see that agreement between the 
two is satisfactory. 

The way to estimate a c in the phase-field simulation is as follows, (i) We consider the 
sharp interface location z — (,{x,t) in the linear analysis to be the null point <p = in 
the phase-field modelling. Numerically, such null points are calculated using the linear 
interpolation of <f> data on the grid points (we use this correspondence throughout this 
paper), (ii) For a given acceleration amplitude a, we monitor the oscillating interface 
position at x = L x /2 as a function of time. We perform simulations by changing a from 
the theoretical a c and estimate the critical value at which the temporal oscillation neither 
decays nor grows. More precisely, we take the absolute relative difference between the two 
peak values of the interface positions around t — 2T V and 4T V . If the difference is smaller 
than 10 -3 , we regard this a as the critical acceleration of the phase-field simulation. The 
largest error shown in table Q] (the case of k = 28.0) would be due to a relatively large 
Ax = L x /N x . 

Here we set the parameters of the phase-field simulation to be e = Az = 1.80 x 10~ 6 
m, 7 = 6.33 x lO^nAg^s, A = 3o-/(2\/2e) (see (p30]t ) and At = 2.5 x 10" 6 for all 
four cases. The initial perturbed interface is written in terms of the phase variable as 

, (z — \L~/2 — bcos(kx)\ 1 .„ „. 

= tanh | L Zl y -^ | , (3.2) 

where b corresponds to the perturbation amplitude of the interface position. Here we 
take b = 3Az = 3e. We recall that in the unperturbed state the Cahn -Hilliard e quatio n 



has the stationary solution <f>{z) — tanh[(z — L 2 /2)/(\/2e)] (see, e.g. Ijacaminl (J199S 
The initial velocity in the phase-field simulation is set to zero. 

Figure [3] shows temporal variations of the interface location obtained by the Floquet 
analysis results (|3.1[) and the phase- field simulation results with critica l a c values. We 



here follow the representation used in lPerinet. Juric fc Tuckermanl ( 20091) . The initial dis 



crepancies between the linear theory and the simulation remain large until one vibrating 
period T v . However, later than that, the phase-field solutions agree satisfactorily with the 
Floquet analy sis result as shown in figure ES In c ontrast, the discrepancies between the 
simulation by IPerinet. Juric fc Tuckermanl (J2009I ) and the linear theory vanish within a 



quarter of the period of the vibration. 

Finally, we comment on the values of the mobility 7 and the thickness of the interface 
e, which are adjustable numerical parameters of the phase-field simulation. The mobility 
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Figure 3. Time variations of the interface position at x — L x /2 of the phase-field simulation 
(circles) and the Floquet analysis result (curves). The interface position is defined as the phase 
null points cj> = 0. The thick arrows indicate the size of e, which is the length scale of the diffuse 
interface, (a) Perturbation wavenumber k — 48.0 mm -1 (subharmonic case); (b) k = 73.0 mm - 
(harmonic case); (c) k — 94.0 mm" 1 (subharmonic case). 



7 is determined as follows. The physical time scales of the various forces in (|2.1[) can 
be dimensionally estimated for a given length scale l as t g = (lo/g) 1 ^ 2 , t v — IqP/t) and 
t a = {Iqp/o) 1 / 2 , which are times scales of the gravity, viscosity and surface tension, 
respectively. The time scale associated with the mobility is in a similar way given as 
t-y = IP by (12.6[) . The mobility value is determined by requiring t 7 < minl^i,,,^} for 
a small length scale, say Iq = e. This means that the interface's relaxation to equilibrium 
occurs faster than other dynamics. Here we have t g = 4.29 x 10" 4 ,t^ = 3.47 x 10~ 5 and 
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t a = 3.34 x 10~ 5 s. To have the condition t 7 < mm{t g , t v ,tcr}, we set 7 = 6.33 x 10~ 8 
kg _1 m 3 s, which yields i 7 = 2.63 x 10~ 5 s. That is, t 7 / min{t g ,t v , t a } = 0.79 here. If this 
ratio is larger than 1.0, our simulation does not agree with the linear theory. Here we 
take the time step to be At = 2.50 x 10 _6 s, which is about an order of magnitude smaller 
than t 1 . The thickness e is determined as follows. We first fix the perturbation amplitude 
to b = 3Az and then do simulations with e/Az = 0.5, 1.0,2.0, ...,5.0. Good agreement 
with the a c of linear theory is obtained for e/Az ^ 3.0. This leads us to choose e = Az. 

4. Simulation beyond the linear regime 

Having validated the phase-field approach in the linear regime, we now move to non- 
linear regimes. Two cases are considered. The first one concerns plume formation, which 
may eventually lead to droplet ejection at long time. However, we do not follow the mo- 
tion until the ejection or pinchoff occurs. The second case is the period tripling state of 
the two-dimensional Faraday waves. 

4.1. Plume formation 
Plume formation of the Faraday interfaces in two dimensions was studied numerically by 



Wright. Yon fe Pozrikidisl ( 20001 ) using the vortex-sheet method, with which we compare 



the phase-field simulation. The presence of a plume implies that the interface profile 
becomes a multivalued function of the horizontal coordinate. Such situations can be 
numerically handled by the phase-field method without ad hoc adjustments. 

Here the comparison is qualitative. There are three reasons for this: the difference 
in Atwood number A = (pb — Pt)/(pb + Pt), the difference in the boundary conditions, 
and with or without viscosity. Regarding the Atwood number difference, the phase-field 
method with the discretization described in section [2~2l is not able to handle high Atwood 
number A. Our simulation works for at most A ~ 0.40 in practice (this issue is discussed in 
the last <J5|). In contrast, the vortex-sheet calculation by I Wright. Yon fe Pozrikidisl ( 20001 ) 



can cope with up to unit Atwood number. Concerning the difference in the boundary 
conditions, there is no rigid boundary in the vertical direction in their vortex-sheet cal- 
culation. However, in our case fluids are contained between two rigid walls. The third 
reason is that the vortex-sheet calculation is inviscid, while the phase-field method is 

viscous. 

In I Wright. Yon fe Pozrikidisl ( 20001 ) . the overturning interface corresponding to for- 



mation of a plume at A = 0.65 and 1.0 is numerically studied. We observe a similar 
overturning behaviour in the phase-field simulation even at as low as A = 0.10, which 
we show in figure |4j The parameters of the phase simulation are: L x = 1.0m, L z — 
1.0m, a = 9.0 x 10" 2 ms~ 2 ,u; = 2.6728 x lO^s" 1 ,^ = 1.0 x 10 3 kgm- 3 ,p t = 8.1818 x 
10 2 kgm~ 3 ,?7t = 0.2Pas,?7b = 0.2Pas,cr = 7.2 x VOr^vaT 1 , g = 0.0, A = 3<r/(2\/2e),7 = 
5.61 x 10 m kg s and e = Az. The force G in (|2.1j) is here G = acosut, oscillating 
around zero. The initial interface is given by p.2[) with b = 0.01m. The time step and the 
number of grid points are At = 1.18 x 10~ 3 and N x x N z — 256 x 256. Here we take the 
same parame ters L T , uj, py,,b,o~ and the same vi bration force G setting (oscillating around 



the origin) as I Wright. Yon fc Pozrikid is (2000). The differences between our simulations 



and theirs are the Atwood number, the vibration amplitude value (ours are nine times 
larger) , the presence of the top and bottom rigid walls and the presence of viscosity. 

The time variation of the interface location at a fixed horizontal point is shown in figure 
|4|a). The interface location first exhibits a rise and then a fall (here at about 1.5T V and 
2.5T V , respectively, where T v is the vibration period); then it grows steeply (at around 
t = 3.0T V ) and reaches a high plateau (between t = 3.4T V and t = 4.4T V ) with a small 



Simulation of 2D Faraday waves with phase-field modelling 11 

oscillation. This variation is similar to that with A = 0.65 o btained with the vortex- 
sheet method (figure 11(a) of I Wright. Yon fc Pozrikidid ( 20001 )). In our simulation, the 



interface location later than t — 5T V (continuation of figur^a) but not shown) goes 
downwards and shows again a plateau with a small oscillation. Then it goes upwards 
and finally comes back to zero around t = 8T V . This later variation canno t be compared 
with the vortex-sheet method result of I Wright. Yon fc Pozrikidisl ( 20001 ) because their 



calculation is stopped before t = 5T V . 

The phase field <f> at t = 3.5T V is shown in figure SJb), which clearly shows that the 
interface turns over and becomes a multivalued function of the horizontal coordinate 
x. Qualitatively similar interfac e structures, called plumes, are shown in figure 11(b) of 



Wright. Yon fc Pozrikidisl (2000)). To check the boundary effect, we also do the phase-field 



simulation with aspect ratio 2 (L x — 2L Z ). The plume shape and its time evolution do 
not change qualitatively, indicating that the effect of the vertical boundary is small on 
the plume formation. This agreement suggests that the phase-field simulation is valid, at 
least qualitatively, in the overturning regime although it cannot handle cases with high 
Atwood number. 

In figure |D^c) and (e), we plot the vorticity field O = d z u — d x w at the plume state. 
Indeed much of the vorticity concentrates inside the interface region. The vorticity along 
the interface in the right half domain has three peaks as seen in the cont our plot of figure 
[He). In contrast, in the simulation of Wrig ht. Yon fc Pozrikidisl (J2000l ). the strength of 
the vortex sheet has at most two peaks in the developed plume state. In the phase-field 
simulation we observe that the vorticity outside of the interface reaches about 40% of the 
maximum vorticity in the interface region. In addition, thin vortex layers are attached 
near the top and bottom walls. In these layers, the vorticity is about 10% of the maximum 
in the interface region. 

We now comment on the choice of the resolution N x x N z — 256 x 256 used here. 
We observe that the squared modulus of the Fourier coefficients \(j>(k,z,t)\ 2 decreases 
exponentially for large k such as \4>(k, z,t)\ 2 ex exp[— 6(z,t)k]. We measure this factor 
S(z,t) and require ma,x Zit [5(z, t)] ^ 2Ax = 2L X /N X . We take the smallest N x in powers 
of 2 satisfying this relation. We then set Az to about the same value of Ax. 

4.2. Period tripling 

The period tripling state of Faraday waves, in which the period of the standing wave 
becomes three times the usual subharmonic period (which is in turn twice the vib ration 
period ), is foun d in the quasi-tw o -dime nsional experiment by iJiang. Perlin fc Schultzl 
( 19981 ) (see also Perlin fc Schulta (20001)') an d in the thre e-dimensional axi s ymm etric 



2000(1 
fitJC 



experiment by I Jiang. Perlin fc Schultzl (|1998l ) and also by iDas fc Hopfingerl (feoOSl ). In 



the two-dimensional vortcx-shect simulation of Wright, Yon & Pozrikidis (2000), period 
tripling is also observed. From their numerical result, I Wright. Yon fc Pozrikidisl ( 20001 ) 
conclude that the period tripling is caused by the nonlinearity of the irrotational flow. We 
here aim at reproducing the period tripling state with the phase-field modelling within 
our feasible Atwood number range A < 0.40, although the above experiments and the 
vortex-sheet simulation are conducted in the case A ~ 1.0. 

By searching the parameter space, which will be discussed later, indeed we observe 
the period tripling state at low Atwood number A = 0.11, as depicted in figure [5] The 
interface shapes at the maxima are shown in figu re EH some of which actu ally differ from 



the previous observations. In the experiment by I Jiang. Perlin fc Schultzl ( 19981 ). the in- 
terface shapes at the three peaks in the tripling state are classified as sharp crest (called 
mode A by them), flat or dimpled crest (mode B) and round crest (mode C), respec- 
tively. The corresponding shapes of the phase- field simulation shown in figure [6^ a-c) are 
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Figure 4. Faraday interface in a nonlinear regime, showing an overturning behaviour in the 
phase-field simulation, (a) Time variation of the interface position at x = L x /2. The interface 
position is defined from the phase field cf> as described in section [3] (b) Grey-scale coded phase 
field <f> at time t — 3.5T V . (c) Grey-scale coded vorticity field at the same instant, (d) Interface 
profiles at four instants, (e) Contours of the vorticity at t = 3.5T V . Contour values are from -20.0 
to 20.0 by 4.0. The solid and dashed lines correspond to non-negative and negative vorticity 
values. 
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qualitatively different except for the round crest (figure Etc)). In particular, the hour- 
glass plume shown in figure Eta) and waistless plume in figure Eth) are different from 
the sharp crest (mode A) and the flat or di mpled crest (mode B). These s harp and flat 
crests are ass ociated with wave break ing in Ijiang. Perlin fe Schultzl (|1998j) but not nec- 



essarily so in Das fc Hopfingerl ([2008). In the phase- field simulation with this parameter 
set, wave breaking does not occ ur. In the period tripling state observed numerically by 
Wright. Yon fc Pozrikidisl ( 20001) . plume- type interfaces are not observed. 



The vorticity contours at the three peaks are shown in figure E(d-f). The number of 
vortex peaks inside the interface region is three for the hourglass plume and one for 
the waistless plume and the round crest. Outside of the interface region the vorticity is 
weaker but not zero. 

To find the parameters of the period tripling state in practice, we use as a guide 
the phase diagram made by Ijiang. Perlin fe Schultzl (|1998l ) in the space of the detuning 



parameter p and forcing parameter q. We search this p, q-space horizontally (constant 
p) by going away from the neutral curve as shown in figure [7] We take seven parameter 
sets along this line. The period tripling s tate is found at the middle point. In the phase 



npln 

tz(l 



diagram in I Jiang. Perlin fe Schultzl (J1998J), the period tripling state extends for the large 



q region. But our case is very different. The parameters of the tripling state are: L x = 
1.46 x 10" 4 m, L z = 2.31 x 10~ 4 m, a/g = 9.5, u = 2ir x 10 2 s _1 , p h = 5.19933 x 10 2 kg 
m~ 3 , p t = 4.15667 x 10 2 kg m" 3 , n h = 3.908 x 10~ 5 Pa s, r) t = 3.124 x 10~ 5 Pa s, a = 
2.181 x 10~ 6 N m -1 , g = 9.8066 m s~ 2 , e = Az = 1.80 x 10" 5 m, 7 = 6.33 x 10~ 8 m 3 kg- 1 s, 
A = 3cr/(2v2e) and N x x N z = 128 x 128. The spatial resolution is determined in the 
same way as described in subsection 14. II 

Now we discuss the behaviours at other parameter points shown in figure [71 For three q 
values smaller than the period tripling value, we observe that the period of the standing 
wave remains twice the vibration period (the subharmonic period). The interface shape 
at the maximum and minimum surface elevations are hourglass-plume type. In other 
words, the temporal symmetry is kept as we approach the period tripling state from the 
left w hile in the experiment and numerical simulation by Ijiang. Ting. Perl in & Sc hultzl 
( 19961 ) the symmetry is broken before reaching the tripling state. This discrepancy may be 



due to the difference in the Atwood number or the difference in the dissipation caused by 
the sidewalls. For the next larger q parameter than that of the tripling state, we observe 
that the period of the standing wave returns to the subharmonic period (2T V ). We do 
not have an explanation for this. The interface shape at the maximum and minimum 
elevation is also of the hourglass-plume type. For the two largest q parameters that we 
calculate, the oscillation becomes close to quasi-periodic. The interface shapes at relative 
maximum elevations in this parameter range take the three forms shown in figure EI a-c) 
almost randomly. 



5. Discussion and concluding remarks 

We have app lied the phase-fi eld modelling of the binary fluids with the Cahn-Hilliard 
equation due to ljacqminl ( 19991 ) to numerical simulations of the Faraday wave problem in 



two spatial dimensions. Here we have solved the Navicr-Stokcs equations for both the top 
and bottom fluids with rigid boundary conditions on the top and bottom walls and with 
periodic boundary condition for the side walls. Validation of this phase-field simulation 
is checked quantitatively in the linear regime and qualitatively in the nonlinear regime. 
In the linear regime, our simu lation agrees quantitatively well with the Floquet analysis 
bv lKumar fc Tuckerman ( 1994) in three differe nt branches, which is the benchmark test 



proposed in lPerinet. Juric fc Tuckermanl (|2009[ ) for the Faraday waves. 
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Figure 5. Temporal variation of the interface position at x — L x /2, showing the period 
tripling. Here T v is the period of the vibration forcing. 
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Figure 6. Interface profiles corresponding to the three peaks of the period tripling state (a) 
t = 19.4T V ; (b) t — 21.4T V ; (c) t = 23.4T V . They correspond to the last three peaks shown 
in figure [5] Vorticity contours (solid and dotted lines correspond to non- negative and negative 
values, respectively) at the same instances: (d) t = 19.4T V , the contour values are from —5.0 to 
5.0 by 1.0; (e) t = 21.4T V , the values are from -7.0 to 7.0 by 1.0; (f) t = 23.4T V , the values are 
the same as (e). 
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Figure 7. The p, q parameter space for seeking the period tripling regime. Here 
p — {2fl(k)/oj} 2 is the detuning parameter and q — 2Aka/co 2 is the forcing parameter, where 
Q 2 (k) = kg[A + ak 2 /{g(p h + p t )}\. The points are q = 2.00, 2.11, 2.17, 2.23, 2.29, 2.35, 2.41 and 
constant p = 2.35. The curve correspo nds to the neutral stability line obtained by the linear 
analysis of iKumar fc Tuckermanl (| 19941 ) (the method described in Sj3]is used). 



In the nonlinear regime, we are not able to validate t he pre sent simulation against 
experimental results unlike IPerinet. Juric fc Tuckermanl ( 20091 ). Nevertheless we have 
considered two cases. Both cases involve plume formation where the two- fluid interface 
becomes a multivalued function of a horizontal coordinate. Such situations can be a good 
testing ground for the phase-field method since the method does not break down when 
the interface turns over. 

The non linear case we considered first i s the bursting plume formation which was 
studied by IWright. Yon fc Pozrikidisl (|2000l ) with the vortex-sheet method. Their sys- 
tem setting (the boundary conditions, the dissipation and close-to- unity Atwood num- 
ber, especially) is very different from the present binary fluid setting. Consequently the 
comparison is necessarily qualitative. The phase-field method with the numerical scheme 
described in section 3 works only in the Atwood number range A < 0.40. Nevertheless we 
found for very low Atwood number A = 0.10 that a similar plu me is formed and the tem- 
poral variation of the interface agrees qualitatively with that in lWright. Yon fc Pozrikidisl 
( 20001 ). 

The second nonlinear case is the period tripling of the tw o-dimensional Fara day waves 

first fo und experimentally bv ljiang. Perlin fc Schultd ( 19981 ) and numerically bv l Wright. Yon fc Pozrikidis 

( 20001 ). Again, in spite of the rather big system differences with these studies, we have 

observed the period tripling state in the phase-field simulation with low Atwood number 

A — 0.11. This result adds further eviden ce for the robustness of the period tripling as 

discussed in ljiang. Perlin fc Schulta (|1998T ). More detailed numerical study of the tripling 

state is currently under way and will be reported elsewhere. It would be interesting to 

perform a laboratory experiment with the same physical parameters as in £14.21 to check 

whether or not the period tripling is observed. 

Concerning the limitation on the Atwood number A < 0.40 in the present phase-field 
simulation, we do not have an explanation for why it is so. For A ^ 0.40, we observe 
that, if we put a larger number of grid points inside the interface region, we can postpone 
the breakdown of the simulation until somewhat later (we can continue the simulation 
for a larger number of time steps). However, increa sing the number of grid points there 
may cause a new problem. Although lJacgminl ( 20001 1 successfully reaches Atwood number 
0.98, he reports that the interface becomes much wider than in reality. It is not clear 
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to us for the present whether or not we can overcome these problems by improving the 
temporal and spatial discretizations of the Cahn-Hilliard and Navier-Stokes equations. 

Finally, we comment on extension of the phase-field method to three-dimensional 
Faraday waves. It is straightforward. We believe that the three-dimensional p hase-field 



metho d, once extended, can be complementary to the simulation method by Perinet. Juric fc Tuckerman 



(2009) since the two methods are different. There is an interesting instance of three- 



dimensional Faraday waves, where the interface becomes multivalued . That is the prop- 
agating solitary state observed by iLioubashevski. Arbell fc Finebera ( 19961 ). Our simu- 



lation code can easily study such a state. 
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